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Abstract. First-order hyperbolic systems are promising as a basis for numerical 
integration of Einstein's equations. In previous work, the lapse and shift have 
typically not been considered part of the hyperbolic system and have been prescribed 
independently. This can be expensive computationally, especially if the prescription 
involves solving elliptic equations. Therefore, including the lapse and shift in the 
hyperbolic system could be advantageous for numerical work. In this paper, two first- 
order symmetrizable hyperbolic systems are presented that include the lapse and shift 
as dynamical fields and have only physical characteristic speeds. 

PACS numbers: 04.20.-q, 04.25.Dm 
1. Introduction 

There has been considerable interest recently in first-order hyperbolic systems for 
Einstein's equations ([||, |^ and references therein). These systems have been used in 
the past to prove that general relativity has a well-posed initial value formulation |^. 
Much of the recent interest is based on the advantages that hyperbolic formulations offer 
to numerical simulations [|], |^ . The main advantage is that imposing physical boundary 
conditions is much easier in the framework of a hyperbolic system than a non-hyperbolic 
one. This is especially true for boundary conditions inside a black hole horizon |^. 
Indeed, if the hyperbohc system has only physical characteristic speeds — that is, if the 
characteristic fields propagate only on the light cones of spacetime or normal to the time 
slices — then the boundary condition inside the horizon on fields propagating into the 
numerical grid has no effect on the dynamics outside the horizon.|| Therefore, in this 
case, any convenient boundary condition can be imposed inside the horizon. This is a 
significant advantage when simulating black holes. 

It is particularly important to come up with stable numerical schemes to evolve 
black holes since simulations of black hole collisions have an important role to play 

J It is sufficient for the characteristic fields to propagate on or within the light cones for this to be 
true. 
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in the detection and analysis of gravitational waves. These simulations will be used 
in several stages of data analysis for gravitational wave detectors such as the Laser 
Interferometer Gravitational- Wave Observatory. First, the simulations are expected to 
yield a bank of gravitational waveforms that will be used to detect the presence of a 
gravitational signal in the detector output. Once a signal has been detected, numerical 
simulations will be used to extract binary parameters such as masses from the signal, 
to test general relativity, and to do other interesting physics. 

Previous numerical work has generally been restricted to systems that do not treat 
the lapse and shift as dynamical fields, but rather take them to be external to the 
system and prescribe them independently. Freedom in choosing these gauge fields 
corresponds to freedom in choosing coordinates for spacetime. This freedom can be 
used for a variety of purposes, e.g., to prevent the occurrence of coordinate singularities 
and reduce coordinate shear and to adapt the coordinate system to the particular 
problem under consideration. When simulating black holes, it is helpful to choose the 
shift so that numerical grid points do not fall into the holes. When simulating binary 
black holes, it may be advantageous to implement gauge conditions which generate 
corotating coordinates ||^, [1^ . 



Some of the favored gauge choices in numerical relativity p, |T^ require solution of 
elliptic equations for the lapse and shift, which is expensive computationally. It would be 
more efficient to evolve the gauge fields as part of the hyperbolic system. However, it is 
important to keep some freedom in choosing the gauge in order to allow the coordinates 
to be adapted to fit specific needs. The purpose of this paper is to present two first-order 
symmetrizable hyperbolic systems which include the lapse and shift as dynamical fields 
and allow four functions of spacetime to be specified freely in the gauge prescription. 

Previous work in this direction includes , in which the authors present a weakly 



hyperbolic system|§] that incorporates the gauge fields in the system, and |T2|, in which 
the authors present a new class of dynamical gauge conditions which are not, however, 
part of a first-order hyberbolic system. 

The first hyperbolic system presented in this paper is based on the work of Fischer 
and Marsden it uses generalized harmonic coordinates and evolves 50 fields. It is 
promising as a basis for numerical work. The second system is based on the work of 
Kidder, Scheel, and Teukolsky and Lindblom and Scheel [0; it evolves 70 fields. 
This system is not practical for numerical implementation. Its main use is theoretical: 
it allows one to show that any solution to Einstein's equations in any gauge can be 
obtained using hyperbolic evolution of the entire metric, including the gauge fields. 
Both systems have only physical characteristic speeds. 

In this paper, Greek indices range over 0, 1, 2, 3 and Latin indices over 1, 2, 3. The 
sign conventions are those of [|T^ with G = c = 1. The analysis of this paper is done 
within the framework of a 3-1-1 split of spacetime (see, e.g., M, O). In this framework. 



§ I refer to the full system including lapse and shift as dynamical fields; if the shift is considered a 
fixed spacetime function and not a dynamical field, then the system becomes strongly hyperbolic. 
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the spacetime metric is expressed as 



and the inverse 4-metric as 




;i.2) 



where a is the lapse, /3* is the shift, is the spatial 3-metric with inverse 7*-', and 
A = TijA- The unit normal to the time slices is denoted by n^. 

I restrict attention in this paper to the vacuum Einstein equations. 

2. System I 

2.1. Fischer-Marsden system 

Let us first briefly review the Fischer-Marsden system for Einstein's equations. 
They employ the 50 fields g^^, k^^ = dtg^u, and di^^ = dig^y. Using harmonic 
coordinates, they reduce the vacuum Einstein equations R^u = to the following first- 
order symmetric hyperbolic system: 



dtg^iu = 

-g'^'dik^, - 2g''%k^,, - g^'ddj^u = 

9''dtdj^,-g'^d,k^, = 0, (2.1) 

where H^^i, is a function of the fields g^^,, k^^, di^^, only and not their derivatives. This 
system is obtained by setting to zero a reduced form of the Ricci tensor that is equal 
to the full Ricci tensor in harmonic coordinates. Using earlier work of Choquet-Bruhat 
( [1T6| and references therein), Fischer and Marsden show that if the initial data for 
satisfy the harmonic coordinate condition and the constraint equations, then the 
solution of ( p.l|) corresponding to these initial data continues to satisfy the harmonic 
coordinate condition off the initial hypersurface. Therefore, a solution of (|2.1| ) is also a 
solution of the vacuum Einstein equations. 

The Fischer-Marsden system ( p.l|) has two drawbacks when considered as a basis 
for numerical integration of Einstein's equations. The first is the restriction to harmonic 
coordinates: this eliminates the freedom to choose coordinates best suited for the 
physical problem at hand. While harmonic coordinates have been used successfully 
in some previous work ([^ and references therein) and are being strongly advocated 



for a wide variety of applications |T^, it has not yet been established whether they are 



useful for simulating black hole collisions, for example. 

The second drawback is that the Fischer-Marsden system has nonphysical 
characteristic speeds. As discussed above, systems with only physical characteristic 
speeds are better suited for numerical relativity, especially for black hole simulations 
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0, 1^. The characteristic speeds of the Fischer-Marsden system can be calculated as 
follows: first write ( p.!]) in the form 

dtu + A\t, x\ u)diU = F{t, x\u), (2.2) 

where u is a column vector composed of the fields {u = {gfj,u, k^^, di^u)"^ for the Fischer- 
Marsden system), and the matrices and column vector F can depend on space and 
time and on the fields but not their derivatives. Pick a unit spatial covector C,i (i.e., 
= 1) ^iid compute the eigenvalues A of the matrix A^^f, A are the characteristic 
speeds in the direction ^j. For physical characteristic speeds, we require A = — 
— ± a (see, e.g., 0). However, the Fischer-Marsden system has A = 0, — ± a. 



2.2. Generalized harmonic coordinates 

In this paper, I modify the Fischer-Marsden system to eliminate nonphysical 
characteristic speeds and generalize it to include a broader range of coordinate systems. 
Let us begin by defining F'^ = (/"^F^q,^ and F^ = g^^T^ , where T'^ ai3 are the Christoffel 
symbols associated with the metric g^j^u and the coordinates x^. The Ricci tensor can 
be written as 



where 



and 



R^y = Rfj_u + V(^Fj,), (2.3) 



R,u = —g'^^d^dpg,, + H^,{g, dg), (2.4) 



H,. = g''^g'"'{d^gj{d(,g,^) - g^'^g'^^T ^p^T.^p. (2.5) 

I generalize harmonic coordinates using Friedrich's gauge source functions H by 
setting 

F'^ = -V„V°x^ = /^(t, x^), (2.6) 

where the coordinates are treated as scalar fields in the expression Vq,V°'x'^, and 
are arbitrary but predetermined functions of space and time. These functions can be 
used to tailor the coordinates to fit specific needs. 

Consider the reduced equations obtained by setting 

+ V(^/,) = 0, (2.7) 

where = g^uf^- Equation ( p.7|) will be used to write down a first-order symmetrizable 



hyperbolic system in section p.3| . Hence we must show that a solution to ( p.7|) yields 
a solution to the vacuum Einstein equations R^u = under appropriate conditions. I 
follow an argument due to Friedrich p, Q which is based on earlier work by Choquet- 
Bruhat (||16| and references therein). 

Let g^i, be a solution to ( |2.7| ). Compute F^ and R^^, from g^^^, and let h^ = T^^ — f^. 
Then R^^, = R^^, + V(^F^) = V(^/ijy) where h^ = V^ — f^. The Einstein tensor is 

R 1 
G^^ = Rfj,u — —giiv = ^{^hi,) — -gfj.v^ah"', (2.8) 
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and the contracted Bianchi identities V^G^'^ = imply 

Vf'Vf.h" + R\h>' = 0, (2.9) 

which is the subsidiary equation derived by Friedrich [|, |]. Since this is a hnear 
homogeneous wave equation for h^, we conclude that if /i'^ = and V^h^ = on 
the initial hypersurface, then = in a neighborhood of the initial hypersurface. 
This imphes Rfj,u = V {^h^) = in this neighborhood. So g^j^u is a solution to the 
vacuum Einstein equations in a neighborhood of the initial hypersurface. This solution 
is obtained in coordinates satisfying = f^. 
We therefore need to ensure 

[F'^ - nt=o = 0, (2.10) 
[V,iT^-n]t=o = 0, (2.11) 

where the time slice t = represents the initial hypersurface. Given a spatial 3-metric 
'jij and an extrinsic curvature Kij that satisfy the constraint equations, we will construct 
initial data for our system such that ( p.lO| ) is satisfied. Equation ( p.ll| ) will then follow 
from the constraint equations. This will be discussed in detail in section [27 



2.3. System I 
Define the fields 

k^^,. = dtg^j^v - P^djgf,^, (2.12) 
dif,u = dig^u- (2.13) 

Here and throughout this section, /3* will be considered convenient shorthand for 
—g'^^/g^'^, and similarly a for (—5'°°)"^''^. The new field k^^j is a replacement for k^j^i, 
and has been introduced to eliminate nonphysical characteristic speeds. 

The first-order symmetrizable hyperbolic system presented in this section is based 
on the 50 fields g^j^u, k^j^^y, and di^^. The definition ( |2.12| ) yields an expression for dtg^u 
in terms of the 50 fields and their first spatial derivatives. An expression for dtdi^^, is 
obtained through equality of mixed partials: dtdi^^ = didtg^u = di{k^y+f3^djfj_„). Finally, 
an expression for dtk^y is obtained from the reduced equation ( |2.7] ). To summarize, we 
have the first-order system 

^tg^,u + -^digf,^ = k^^, (2.14) 



dtk^u H -^dik^y H —didj^^ = —g "di^i^k^j H —[Hf^i, + <9(^/;^) — r°^!//a]; (2-15) 



dtdi^ii ~\- -^^djdi^y dik^p g dj^yd^^]^.^ (2.16) 

where 7*-' = {g^^)~'^{g^^g^-' — g^^g^-') is the inverse of the 3-metric jij. In ( |2.15|) , if^j. 



is to be expressed via ( p.5|) in terms of the fields only and not their derivatives [using 
( p.l2| ) and ( p.l3|) ]. In addition, in (|2.14|) -( ^.16|) , the inverse 4-metric is considered to 
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be a function of (yf^i, and not a fundamental field. In deriving these expressions, I have 
used the relation 



dag'"' = -g^'g^'^dagex. 
The system ( |2.14| )-( ]2.16| ) will be called system I. 



(2.17) 



2.4- Initial data 



It remains to specify how to set initial data for system I to ensure ( |2.10 ) and ( 2.11| ) 
are satisfied. Begin with a solution (7jj, Kij) of the constraint equations, where Kij 
represents the extrinsic curvature of the initial hypersurface. First set gij = 'jij. We are 
free to choose go^^ on the initial hypersurface as long as goo < goigojY'' ■ This requirement 
is equivalent to > and implies (7°° < 0. Freedom in choosing gof_i corresponds to 
freedom in choosing the lapse and shift at t = 0. 

We now have g^u\t=o- Next set c/j^,^ = dig^^\t=o- The final step is to fill in k^^^, from 
Kij and the requirement ( |2.10| ). The extrinsic curvature can be expressed as 



K., 



From this we deduce 



(2.18) 



(2.19) 



kij = -2aKij + 2gk{idj)l3^, 

which can be used to fill in kij at t = 0. 

The quantities fco^ are obtained from the requirement ( p.lO| ). Writing out in 
terms of the metric and its first derivatives, we obtain 



(2.20) 
v^■^(2.21) 



F° = -a-^{dta - (3'dia + a^K), 

r = -a-\dtf3' - (3^djf3') + a-%dta - p^dja + a^K)p' - a-^^d^a + ^^^V^kT 

where K = Y^Kij, and *^^'*F*jfc are the Christoffel symbols associated with the 3-metric 
■jij and the spatial coordinates . Setting F^ = gives us expressions for dta and dtP^ 
which we use to fill in /cq^ at t = 0: 



where 



koi 
Bi -- 



Bi + (3'' kij, 

: 2a%af + K) + 2^'Bi + ^'^^k 



-a\gi^r + a-^dia-^^^Vijkl- 



(2.22) 
(2.23) 

(2.24) 



The initial data for system I is now complete and satisfies the constraint equations 



\t=o 







(2.25) 



and the requirement (|2.1CI|) . This in fact implies that the requirement ( |2.11|) is satisfied. 



The argument follows earlier work |T6[ on the reduction of Einstein's equations using 
harmonic coordinates. From (|2.8| ) and (|2.25| ), we deduce 

2rfVi^^K) - Uf.Vah'^ = 0. (2.26) 
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Here and in the remainder of the paragraph, all quantities are evaluated at t = 0. We 
know /i^ = — = on the initial hypersurface, so v^'Wyh^ = for any spatial vector 
v'^ (i.e., for v'^ satisfying f'^n^ = 0). It remains to show n^Vph^ = 0. By contracting 
( p.26| ) with v^, we obtain v^n^Vuh^ = 0. Furthermore, Vah'^ 
Contracting ( p.26[ ) with n^, we obtain n^n'^'W ^hy 
and so ( ^.111 ) is satisfied. 

Therefore, a solution {g^u,k^y,di^y) to system I with initial data as constructed 
above yields a solution g^j^y to the vacuum Einstein equations. 



0. It follows that n'^Vyh'' = 



2. 5. Hyperholicity of system I 

System I is symmetrizable hyperbolic. To see this, let u 
and write equations ( p.l4|) -( ^.16|) in the form 
A' to be 



d\^Pi d2^pi c^s/xi^) 



121). This determines the 50 x 50 matrices 



A' 



f 


-(3'I 
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-82'! 












V 















J 



(2.27) 



Here and in equations (|2.28| )- (|2.29|) , is the 10 x 10 zero matrix and / is the 10 x 10 
identity matrix. It can be checked easily that the positive definite symmetric 50 x 50 
matrix 



H 



/ / 














\ 
























7"/ 


712/ 


^13/ 










7I2/ 


7^2/ 


723/ 




Vo 





^13/ 


7^3/ 


733/ 


/ 



(2.28) 



is a symmetrizer for the system, i.e., HA^ are symmetric matrices. 

Moreover, system I has only physical characteristic speeds; that is, the eigenvalues 



of A^i are A_ 



-p'^i + a, \o = -P'^i, and A_ 



'P^^i ~ Let rii and Xi be unit 



spatial covectors that form an orthonormal triad with ^j, that is, rf^i = = = ^^Xi 
and ri^rii = 1 = where 77* = 7*-^77j and = Y'^Xj- Let us construct a 50 x 50 matrix 
V whose columns are 50 linearly independent eigenvectors of One such matrix is 



V 





—al 

V 6/ 



/ 











r]il 
V3I 



\ 



(2.29) 




al 
Xil U 

X2I 

X3I 6^ J 

The first ten columns of V are eigenvectors of A^C,i with eigenvalue A+; the next 
thirty columns have eigenvalue Aq; and the last ten columns have eigenvalue A_. 
The characteristic fields in the direction are obtained from V~^u and are given by 



kfiu i ^ difj^y^ 9fiui V difj^y^ and x ci; 
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3. System II 

In this section, all indices are lowered and raised by the spatial 3-metric jij and its 
inverse 7*-^. The second system presented in this paper is based on a hyperbolic system 
in 0, which is in turn based on the ADM equations [|19|. The system in 0, called 
system 1, employs the 30 fields ■jij, Kij, and 

dkij = dk'jij. (3.1) 

It is obtained by densitizing the lapse and adding multiples of the constraint equations 
to the evolution equations. The relevant constraints are the Hamiltonian constraint 

C =^{^^^R~KijK'^ + K^) = 0, (3.2) 
the momentum constraints 

Ci = DjKi^ - DiK = 0, (3.3) 
and the constraint 

Cijki = d[idj]ki = 0, (3.4) 
where ^^^R and Di are the Ricci scalar and covariant derivative associated with 7jj, and 

System 1 has five free parameters that govern how to densitize the lapse and how 
much of the constraints to add; these parameters determine the system's hyperbolicity. 
In fact, it has been shown [T^ that for a certain range of these parameters, system 1 is 



symmetrizable hyperbolic and has only physical characteristic speeds. 

Here I construct a first-order symmetrizable hyperbolic system based on system 1 
that includes the lapse and shift in the system. Let us begin by defining the densitized 
lapse 

Q = ln(a7-'/'), (3.5) 
where 7 = det(7jj). Next define the new fields 

Qi ^iQ : Qij diOjQ^ 

W = diP^, 6,/ = d,d,(3'. (3.6) 

Note that Qij = Q{ij) and bij'' = b(^ij)^. The hyperbolic system presented in this section 
is based on the 70 fields 7,^-, Kij, dkij, Q,Qi, Qij, l3\ bi^ , hj^ . 

Expressions for time derivatives of these fields are obtained as follows. First, dt'^ij 
is obtained from (|2.18|) : 

dtli^ - (3^dkiij = -2aKij + 27^(^6^-)^ (3.7) 

This is one of the ADM evolution equations with the new fields (|3.6| ) substituted in. 
Here and henceforth, it is understood that a is to be rewritten in terms of Q using (|3.5|) . 
Following 0, I add dajijC and C2tt7'""Cm(jj)n to the second ADM evolution equation 
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(which is equation (2.9) in 0), where Ci and (2 are free parameters. Rewriting this 
equation in terms of the new fields (|3.6|), we obtain 



dtKij = (3'^dkKij — -a^'^'^[dyndnij + '2'd(^idj)mn — (1 — C2)d(id\mn\j) — (1 + C2)dmd(ij)n 

- Cilijl^Kdmdkin - dkdimn)] + '^Kk{ibj)'' - a[2Ki^K'^ j - KKij + Q^j 

1 ~ S 

_l_ / J J \( d™" d™" 4- rl^^ d —d d _1_ 

~r \f-l'(^ij")m 2^'-'"mij J y'-'' ^; J T w i'-''[nm]j ^^imn'-''j ~r ^i^j 

+ QcjC^i) + -^didj] + -C,iOilij{dmd"^ — dmd"^ — -^dmd^ — -dkimd"^^^ 

+ -^dki^d^'"^ - K^nK"''' + K^), (3.8) 

where (ij = ■j^'^dijk and (ij = 'j^^djki. 

Using equahty of mixed partials, we have dtdkij = dkdt'jij which, together with a 
spatial derivative of (p.7|), yields an evolution equation for dkij- Following |0, I add 
C3Q^7fc(jC*j) and Ci^^lifik to this equation and use (|3.6| ) to obtain 

dtdkij = P'^dmdkij + a7'""[C3(7fc(i<9|mif„|j) - 'Jk{idj)Kmn) + Ulij{dmKnk - dkK^n)] 

~l~ dfjiijbk ~l~ 2dkm(ibj) 

- aK,,{2Qk + dk) + aUnAKk^i^d"^ ' dT) + ^^"^"4™™] 

+ aCz[lk(^K,)^{U^ - dJ^) + \K^''lki^d,)mnl (3.9) 

where Cs and ^4 are free parameters. The parameters (Ci, C25 Cs; C4) i^i the above equations 
correspond to the parameters (7,C)''?)X) i^i II]- The parameter a in has been set to 
1/2 by the definition (|3j). 

The next step is to specify evolution equations for the lapse density and shift. 
Spatial derivatives of these equations will then yield evolution equations for the fields 
( ^161) . I consider a particular form for the lapse density and shift evolution equations, a 
form that results in a symmetrizable hyperbolic system but yet allows four functions of 
spacetime to be freely specified. The equations are 

dtQ- P'diQ = ^\t,x^;Q), (3.10) 
dtf3'- 13=0,(3' = #(t,x^Q,/3'"), (3.11) 

where ip^ are arbitrary but predetermined functions of space, time, and lapse density 
(and of shift in the case of ip^). 

Evolution equations for the fields ( |3.6| ) are obtained by taking spatial derivatives 
of ( ^.1(J| ) and ( |3.11| ), and using equality of mixed partials. For example, dtQi = didtQ = 
di{(3^djQ + We obtain 

dtQi - (3'djQi = QjW + d,i;\ (3.12) 
dtQi, - (3''dkQij = 2Q,(i6,)^ + Q^fe,/ + didj^p^, (3.13) 
dt(3,= - (3'dk(3,= = b,%= + d,^=, (3.14) 

dtb,," - rdjo^i" = 2b^rh)j + ^^r^m" + did,ip\ (3.15) 
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where it is understood that the spatial derivatives of ip'^ are to be written, using ( |3.6[ 
in terms of fields only and not derivatives of fields. 



When the system ( p.7|) -( pTT5D , called system II, is put in the form ( |2.2| ) with 

u = {jij, Kij,di:ij,Q,Qi,Qij, P\bi^ ,bij^)'^ , the 70 x 70 matrices have the block 
diagonal form 

^.^M'SOXSO 030X40 \ (3^^g^ 
V ^40x30 ~H -'40x40 / 

The nontrivial parts A* of A* come from the evolution equations ( |3.7| )-( |0| ) for the 
30 fields •jij, Kij, dkij. Since the principal parts of these equations are identical (after 
relabeling the free parameters as indicated above) to the principal parts of the system 
1 evolution equations for •jij, Kij,dkij given in the matrices A'- are identical to the 
corresponding matrices in 0. This implies that if system 1 is symmetrizable, so is 
system II. Indeed, the matrix 

}J — ( -^30x30 030x40 j ^^-^ 
\ 040x30 -^40x40 / 

where if30x3o symmetrizes system 1, is a symmetrizer for system II. In other words, if 
the 30 X 30 matrices HA^ are symmetric, then so are the 70 x 70 matrices HA\ In 
addition, the characteristic fields Q, Qi, Qij, P\ bi\ bij^ all propagate normal to the time 
slices. 



It has been shown |]13[ that system 1 in |P] is symmetrizable and has only physical 
characteristic speeds when the free parameters are chosen as follows: 

^8 _ 4 + lOCi + 4C2 + 6C1C2 

5 + IOC1 + 7C2 + 6C1C2' 5 + IOC1 + 7C2 + 6C1C2' 

-5/3<C2<0, 5 + lOCi + 7C2 + 6C1C2 ^ 0. (3.18) 

We conclude that for the same choice of parameters, system II is symmetrizable and 
has only physical characteristic speeds. 

System II is not practical for numerical implementation. Since the lapse density and 
shift evolution equations (|3.10|) and ( p.ll|) decouple from the rest of the system, they can 
be evolved separately to obtain the lapse density and shift as spacetime functions. These 
functions can then be substituted into system 1 in 0|. Therefore, the full seventy-field 
system II does not need to be evolved; the thirty-field system 1 suffices. 

However, system II is useful from a theoretical point of view. Consider a solution 
of Einstein's equations in an arbitrary gauge. Using the densitized lapse and shift 
from this solution, compute the left-hand sides of equations ( |3.10| ) and p.ll| ). Set the 
spacetime functions ip^ equal to these computed quantities. Take initial values for the 
fields in system II from the spacetime metric under consideration. System II can now 
be used, with these initial values and with as defined above, to obtain the entire 
metric by evolving hyperbolic equations that are part of a symmetrizable system with 
only physical characteristic speeds. So system II can be used to obtain any solution 
of Einstein's equations in any gauge using hyperbolic evolution for the entire metric. 
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including the densitized lapse and shift. Note, however, that the lapse is not evolved 
directly in this system; it is obtained from the densitized lapse via equation ( |3.5| ). 

4. Future directions 

An important future research direction is to study and understand the stability of 
numerical implementations of system I. It has been shown in previous work that some 
hyperbolic systems are more stable than others when used to simulate black holes in 
three spatial dimensions. The reasons for this behavior are not yet understood. Another 
future research direction is to explore how to use the free functions in system I to 
control the coordinate system. 
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